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Numerical simulations indicate that the Born rule does not need to be postulated in the de 
Broglie-Bohm pilot- wave theory, but arises dynamically (relaxation to quantum equilibrium). These 
simulations were done for a particle in a two-dimensional box whose wave-function obeys the non- 
relativistic Schrodinger equation and is therefore scalar. The chaotic nature of the de Broglie-Bohm 
trajectories, thanks to the nodes of the wave-function which act as vortices, is crucial for a fast 
relaxation to quantum equilibrium. For spinors, we typically do not expect any node. However, in 
the case of the Dirac equation, the de Broglie-Bohm velocity field has vorticity even in the absence 
of nodes. This observation raises the question of the origin of relaxation to quantum equilibrium for 
fermions. In this article, we provide numerical evidence to show that Dirac particles also undergo 
relaxation, by simulating the evolution of various non-equilibrium distributions for two-dimensional 
systems (the 2D Dirac oscillator and the Dirac particle in a spherical 2D box). 



I. INTRODUCTION 

In the de Broglie-Bohm pilot- wave theory [1-3 , each 
element of an ensemble (consisting of a single particle) 
is described by a wave-function ij^it^x) (i?(t, x)e*'^*^^'^^/^ 
in the polar representation) but also by a position 
x(t). The wave- function always evolves according to the 
Schrodinger equation 

ihdt^pit, x) = -h^m-^Aipit, x) + V{x)ip{t, x) , (1) 

whereas the evolution of x{t) is determined by the stan- 
dard guidance equation 

v{t)=m-'VS{t,x)\^^^^^^. (2) 

The de Broglie-Bohm theory reproduces the predictions 
of standard quantum mechanics provided that the parti- 
cle positions are distributed according to 

p{t,x) = mt,x)f (3) 

over the ensemble. De Broglie's dynamics (defined by Eq. 
([T]) and Eq. ([2|, by contrast to Bohm's dynamics [4 ) en- 
sures that this condition will hold at all times, if it holds 
for some initial time to, a property which is referred to as 
quantum equilibrium [5 , 6 or equivariance [7 . There are 
(at least) four different attitudes towards this condition 
(and some of them are not mutually exclusive): Valentini 
argues that it arises dynamically [5]|6l|8], Diirr, Goldstein 
and Zanghi argue that we happen to live in a typical uni- 
verse [7] , Bricmont takes it as a third postulate [9 , while 
Wiseman suggests that it can be derived as a Bayesian 
prior from a principle of indifference [10]. In this arti- 
cle, we will be concerned with the first attitude, which 
opens a door to a possible new physics, that of quantum 
non-equilibrium (see for instance [11^ and [12 ). 
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If this view is considered seriously, one has to explain 
why we don't see quantum non-equilibrium today. This 
is done by invoking the idea of relaxation to quantum 
equilibrium: if some non-equilibrium distributions ex- 
isted in the past, they are quickly driven dynamically 
to quantum equilibrium. Various numerical simulations 
for two-dimensional systems ([T3HT5]) have shown that 
non-equilibrium distributions rapidly relax to quantum 
equilibrium on a coarse-grained level, provided that the 
wave-function has enough complexity. By complexity, 
we mean that two trajectories originating from neigh- 
borhood points should quickly diverge (which is referred 
to as chaos in a loose sense). The existence of nodes has 
been first connected to chaos in [16]. So when we say 
that the wave-function needs to be complex enough, we 
mean that one needs to superpose a few energy modes 
in order to get a few nodes. Nodes are also associated 
to vorticity: the circulation of the velocity field around 
a closed curved can only be non-zero if there is a node 
inside the curve. More recently, a better understanding 
of chaos, through 'nodal-point-X-point complex' (where 
the X-point is an unstable point associated to the node), 
has also been gained pT] , 

The relaxation simulations have been performed for 
non-relativist ic scalar particles. No simulation has ever 
been done for fermions. And there is a good reason to 
do simulations for fermions, because we typically don't 
expect any node for spinors, no matter how many energy 
modes are being superposed. The reason is the following. 
For scalar particles, we have nodes where 9lt{ip) = and 
3m{tlj) = 0. In 2D, nodes are located at the intersection 
of the curves defined by the two previous equations and 
we typically have a few nodal points (in 3D, we would 
have nodal lines). For a positive-energy spin-^ fermion, 

described by a Pauli spinor ^^^^ ^ nodes are located at 

the intersection of four curves (or surfaces): 9^e(V^i) = 0, 
9^e(V^2) 0, am(V^i) and am(V^2) 0. Therefore 
we typically don't have any node in two or three dimen- 
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sions. This raises the question of relaxation to quantum 
equilibrium. At first, it seems problematic for spinors. 

But the velocity field for spinors (for Dirac particles 
for instance) is essentially different from the standard 
velocity-field for non-relativistic scalar particles: it has 
vorticity even in the absence of nodes. If vorticity is the 
crucial ingredient for relaxation to quantum equilibrium, 
we expect relaxation to quantum equilibrium for Dirac 
particles too. Actually a previous work by Colin and 
Struyve |14 gives some support to that idea. There the 
evolution of non-equilibrium distribution in a 2D square 
box is simulated, the wave-function being nodeless (ex- 
cept at the boundaries). Relaxation is very poor in that 
case, except if 'artificial vortices' are added, such that 
the alternative (and non-standard) velocity field exhibits 
vorticity even in the absence of nodes. This is always pos- 
sible and is referred to as the non-uniqueness problem of 
the de Broglie-Bohm theory [10l[T8]. This example can be 
a toy-model for spinors; however one can argue with rea- 
son that these non-standard velocity fields are unnatural 
and that this example lacks features which are inherent 
to spinors (for instance the spinorial density can't vanish 
at the boundaries). 

That is the main question that we will address in this 
paper. In order to do that, we will consider 2 two- 
dimensional systems for 'confined' Dirac particles: the 
first one is the so-called Dirac oscillator, the second one 
is a Dirac particle confined in a spherical box. In each 
case, we study a superposition of positive-energy solu- 
tions. We will show that the first system, in which 8 
modes are superposed, exhibit fast relaxation, while the 
second system, in which 6 modes are superposed, also 
undergoes relaxation to quantum equilibrium (although 
it is a slower one) . In the second system, we also give an 
example of a non-trivial spinor for which there will never 
be complete relaxation. 

This article is organized in the following way. Section 
II contains the basic details about the Dirac equation in 
3 + 1 and 2 + 1 dimensional spacetimes, as well as the 
corresponding pilot-wave theory. In section III and IV 
we study relaxation to quantum equilibrium for Dirac 
spinors. We conclude in Section V. 

We use units in which h = c = 1. 



II. PILOT- WAVE THEORY FOR A DIRAC 
FERMION 

In this section, we recall the essential properties of the 
Dirac equation in spacetimes of dimension 4 and 3, and 
the corresponding pilot- wave theory. 



A. 3 + 1-dimensional spacetime 



The Dirac equation is 

d ^ ^ - ^ 

i — ilj{t^ x) = —id ■ Vilj{t^ x) + mPip{t^ x) , (4) 



where the matrices aj and ^ must be Hermitian and 
must satisfy the relations {aj^ak} = 25 jk^ = 1 and 
{a„/3} = 0. 

The Dirac equation can be rewritten in a covariant 
form by introducing the matrices 7^ = /3 and 7-^ = j3aj'. 



{iYd^-m)^{t,x) = {) . 



(5) 



The 7-matrices satisfy the relations {7^,7^} = 277^^, 
where r]^^ = diag(l, —1, —1, —1) (Clifford algebra). The 
Pauli-Dirac representation 
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1 
-1 



,7 



di 

(Ti 



and the Weyl (or chiral) representation 
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1 

1 



ai 



(6) 



(7) 



are two common choices of representation for the 7- 
matrices. In this article we use the Weyl representation. 
The conserved 4-current is given by 



(8) 



where ip = ip^j^. 

The positive and negative-energy plane- wave solutions 
are denoted by u{p)e-'^^P^^^'P-^ and ^(p)e^^^^)^+^^-^, 
where E{p) = Y^^^P + rr?. In the Weyl representation, 
if one introduces the right-handed and left-handed eigen- 
states of helicity Xr{p) Xl{p)^ the plane- wave solu- 
tions can be given by 



(9) 
(10) 

(11) 
(12) 




Dirac spinor as ip 



In the limit m 



0, 



There is an easy way to remember this. We write the 

the Dirac equation reduces to a pair of Weyl equations 
for i/jR and i/jl- The Weyl equation for ipR admits 
right-handed positive-energy solutions and left-handed 
negative-energy solutions, whereas the Weyl equation 
for ipL admits left-handed positive-energy solutions and 
right-handed negative-energy solutions. That is indeed 
what we recover from the previous equations in the limit 
E ^p. 
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In the corresponding pilot- wave theory (first proposed 
by Bohm in a reply to Takabayashi [19 ), the particle is 
described by a 4-spinor together with its position 

x{t). The position of the Dirac particle evolves according 
to the guidance equation 



v{t) 



Jit, x) 



. (13) 



Again the choice of the guidance equation ensures that 
the relation 



p{t,x) = \tlJa{t, 



(14) 



a=l 



will hold for any time t provided that it holds for some 
initial time. 

Any positive-energy plane-wave solution of momentum 
p , whether right-handed or left-handed, moves with ve- 
locity p/E whereas negative-energy plane- wave solutions 
move in the opposite direction of momentum. 

Fig. [l] shows the kind of trajectories predicted by the 
theory. The guiding Dirac spinor in this case is a super- 
position of three positive-energy right-handed plane- wave 
solutions 

^x) = 4^(^ii^(pi)e-^^^*+^^^-^>e^^ii^(p2)e-'^^*+'^"'-^^ 



e'^UR{ps)e- 



-iEst-\-ip3-x 



),(15) 



where 



Pi = (1, 0, 1), P2 = (-1, -2, -1) and P3 = (1, -1, 1) 

(16) 

and ur is defined at Eq. (|9|. 



B. 2 + 1-dimensional spacetime 

In a 2 + 1-dimensional spacetime, the Dirac equation 



d d d 

i—^l){t,x,y) = {-iai—-ia2—^mp)ilj{t,x,y) . (17) 

We can take ai = <Ji, 0^2 = (J2 and /3 = (73. 
The covariant form is 



where the 7-matrices are defined by 

7° = (73, 7^ = (73(71, 7^ (7c 

For the plane- wave solutions, we have 



u{Px,Py) 

and 

v{Px^Py) = 



2E 



1 

E-\-m 



-iEt-\-ipxX-\-ipyy 



(18) 
(19) 

(20) 



E -\- m /zPE±iPy.\ 



Trajectories of Dirac particles 
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FIG. 1. Trajectories of Dirac particles in a 3 + 1-dimensional 
spacetime. The guiding spinors are obtained from the one 
given in Eq. (15) by using different values of the mass (3, 6 
and 9). The particles all start from the origin and t G [0, 200]. 



which are respectively positive and negative-energy solu- 
tions. 

In Fig. |2j we plot a few trajectories of particles guided 
by the spinor 

V^(t,f) = ^(ii(gi)e-^^^*+^^"'^-^ + e^'^2i(g2)e"^^^'+^^"'^^ 
v3 

+e^^ii(^3)e"^^^'+^^"^^-^) , (22) 

where 

= (1, 0), 92 = (-1, -2) and 93 = (1, -1) , (23) 

for different values of the mass. 

The circulation of the velocity field 



V ■ dl 



(24) 



can be non-zero, indicating vorticity, even in the absence 
of nodes. Indeed the non-relativistic limit of the Dirac 
equation gives a correction (spin-term) to the standard 
de Broglie-Bohm velocity field [20]. 



III. RELAXATION TO QUANTUM 
EQUILIBRIUM FOR THE 2D DIRAC 
OSCILLATOR 



The equation for the Dirac oscillator is obtained by 
substituting^ by p— zma;/3r in the free Dirac Hamiltonian 



H = a ■ {p — imujpf) + m/3 . 



(25) 



Trajectories for Dirac particles in a 2 + 1-dimensional spacetime 




FIG. 2. Trajectories of Dirac particles in a 2 + 1-dimensional 
spacetime. The gui ding spinors are obtained from the Dirac 
spinor given in Eq. ([22]) by using different values of the mass 
(3, 6 and 9). The particles all start from the origin and t G 
[0,200]. 



This system reduces to a standard oscillator plus a strong 
spin-orbit coupling term in the non-relativistic limit and 
some of its generalizations [22l [23] are relevant to the 
study of quarks confinement. 
The 2D case, with Hamiltonian 



^ ^ I Px-iPy 

\^Px + iPy 



X — iy^ 
X — iy 



was treated by Villalba in [24]. We have verified Vil- 
lalba's construction and we arrive at the same energy 
eigenstates expect that a few intermediate equations are 
different, which is why we include the derivation here. 
Also, it will be useful for the case of the Dirac particle in 
a 2D spherical box. 



A. Energy eigenstates 



We first multiply Eq. ( 26 ) by z and we note that 



dx + ^dy = {dxO + idyO)de + {d^r + idyr)dr 
= r-\-sm0^icosO)de^e'^dr 

= ir-^e'^de^e'^dr , (27) 
dx - idy = {dxO - idyO)de + {d^r - idyr)dr 
= r-^{-smO -icosO)d0 + e-'^dr 

, + e-'^dr . (28) 



-ir 



Therefore we have that: 




iH = 



—muo 



+ e'^dr 

re" 

-re' 









+ (2.) 



The previous operator is applied on 



If we intro- 



duce 



1 _-e 1 -e 

ipi -^e~'^ilj[ and -^e'^i/j^ , 



(30) 



we have that 



lie 



{ir e 



1 



{—ir e 



(31) 



(32) 



Now we look for energy eigenstates: we make the replace- 
ment 



EtikO 



ilj'{t,r,0) ^ilj'{r)e-'^^e 



(33) 



where k is half-integer. From Eqs (29-33), we find that 



ke 

i{E — m)il)\ = (9r H muor)il)2 , 

r 

ka 

i{E + mWn = {dr-^^ mujrU[ , (34) 
r 

from which the following equation for can be deduced: 
-{dl + 4(1 - fc) + mLo{l + 2k) - mWr^)i;[ . (35) 



Now we perform the following change of variables 



= mujr , 



(36) 



which implies that 



d^ilj[ = dr{2mujrd^ilj[) = (2mcja^ + 4m^cjV^a|)V^; • 



(37) 



Eq. (35) becomes 



-{20^ + 4^a| + -(1 - fc) + (1 + 2fc) - Oi^'i , (38) 



where 



muo 



If we start from the ansatz 

= exp{-^/2)CV{0 , 



(39) 



(40) 
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Eq. (38) reduces to an equation for generalized Laguerre phases given in the fohowing table: 



polynomials (7^(0 ^ ^^^(0) 



(41) 



provided that a = k/2 (or a = 1/2 - k/2), fi = k-1/2 
(or -A; + 1/2) and 



4n (or = 4n - 4A: + 2) . 



(42) 



The choice of a (which in turn fixes /i and A^) is fixed 
by the sign of /c, because the wave- function needs to be 
regular at the origin. 

We can always choose units in which h = 1^ c = 1 and 
m = 1. The frequency u remains a free parameter. The 
following table contains the eigenstates of lowest energy 
for the case uj = 1. 



n = 1 
k = 1/2 



n = 1 
k = -1/2 



n = 2 
k = 1/2 



n = 2 
k = -1/2 



^1 



V^\/5-\/5 
2^^/5-y5(l+^/5) 



^1 = 



r(2-r^) 



'/^i(r^-l) 



^2 



^(r^-4r^+2) 



V2 



zr(-r^+2) 



V^2 = 



V2^Vl3-Vl3 
^^*e-"'/"6i(-r^+4r^-2) 
2^^/l3-yT3(l+^/l3) 



n = 1 
A: = 3/2 



n = 1 
A: = -3/2 



n = 2 
A; = 3/2 



n = 2 
A; = -3/2 



V^2 



2e 



7FV5-V5(l+^/5) 



V^2 



r"(-r"+3) 



-^e-^^*e-^'/^6ir(r^-2) 



Try 13-\/l3(l+\/l3) 



e e r(r — 6r +6) 



3V27r 



3V27V 



^1 
^2 



8r^ + 12) 



iVl7t 



8ir(-r^+6r^-6) 



/e^y/ 17- VT7(l+yi7) 



(43) 



(44) 





JT licLot; 




^i4.869 


[L, -l/Z) 


^il.049 


(2 1/2) 


gi4.291 


(2,-1/2) 


gi3.066 


(1,3/2) 


giO.188 


(1,-3/2) 


gil.288 


(2,3/2) 


giO.219 


(2,-3/2) 


gi4.706 



(45) 



We consider 5 different non-equilibrium distributions at 
time t = 0. The first one, denoted by po{t = 0,r, ^), is 
defined by 



2^cos2(^) 



with i^o = 4 , (46) 



if r < Rq^ and zero otherwise. The four remaining distri- 
butions (pi,p2,P3 and P4) are built from po, except that 
they are respectively centered at xi = (2,0), X2 = (0,2), 
xs = (—2, 0) and X4 = (0, —2) and that is substituted 
by = 2. Explicitly we have that: 



2^cos^(^) 

i?2(7r2-4) 



Pj{t = 0, x) 



with = 2 , (47) 



if di{x) = y/{xj{l) - x(l))2 + {xj{2) - x(2))2 < i?, and 
zero otherwise. 

The evolution of the non-equilibrium distribution can 
be computed using a method proposed by Valentini and 
Westman [13 and also used later in [14 and [15 . This 
method relies on the fact that the following ratio 



p{t,x) 



^^t,x)^{t,x) 



(48) 



is preserved along a trajectory. Therefore, if we want 
compute p{t,x)^ we first 'backtrack' the initial position 
(that is, we find the initial condition (^0,^0) which gives 
(t,x) as final position). Then we use Eq. (48) in order 
to find that 



p{t,x) ^l)\t,x)^l){t,x) 



p{to,Xo) 



'0t(to,fo)V^(^o,^o) 



(49) 



For this simulation, we have used a lattice of 2048 x 
2048 points, uniformly distributed inside the box [—5, 5] x 
[—5,5], of coordinates 



B. Simulations 

For the wave-function, we consider a superposition of 
the 8 previous energy eigenstates with equal weights and 



X = -5 + 10/4096 + 10A:/2048 
^ = -5 + 10/4096 + 10//2048 , 



(50) 



with k,l e {0, 1, 2, ... , 2047}). Note that ah the lattice 
points are inside the box (none of them lies on the bound- 
ary of the box), and that is why 10/4098 appears in the 
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expression (and not 10/2048). Each lattice point (repre- 
senting the final position of a particle) is backtracked to 
its initial position by solving the differential equation 



dx{t) ip^ (t, x)dip{t, x) 



dt 



^^t,x)^{t,x) 



(51) 



numerically using the Runge-Kutta-Fehlberg method 
[2F. We have imposed a precision on 0.001 on the back- 
tracked positions. Lattice points that give an insufficient 
precision, or can't be backtracked before a given maxi- 
mum number of iterations (for instance 100000), are dis- 
carded. Hereafter we refer to these lattice points as bad 
ones (or good ones if the converse is true). 

The next step is the coarse-graining, which is done 
by dividing the box [—5, 5] x [—5, 5] in a certain num- 
ber of non-overlapping coarse-graining cells, and averag- 
ing the densities corresponding to the lattice-points con- 
tained in each coarse-graining cell. If we divide the box 
in 32 X 32 non-overlapping coarse-graining cells, the mean 
percentage of good lattice points per coarse-graining cell 
is 98.47% and the worst ceh has 85.06%. 

For the plots, we have actually used a smooth coarse- 
graining of the densities (the same that was used in [14 ) 
which is similar to a standard coarse-graining except that 
it is defined over overlapping CG cells. In this case, all 
the overlapping CG cells can be obtained from a square 
cell of side 10/16 located in the lower left corner by trans- 



lating it by m 



10 



in the up-direction, and by n 



10 



16x8 "^^^ c.x^ "^16x8 

the right-direction, where m^n G {0,120}. The smooth 
coarse-graining is denoted by 





p{t,x) or (V^t^)(t,x) 



(52) 



FIG. 3. Density plots for t G {0,50,100}. The first column 
displays the evolution of the smoothly coarse-grained stan- 
dard density (ip^ip), while the second column displays the evo- 
lution of the smoothly coarse-grained non-equilibrium density 
(po)- The smooth c oars e-graining is defined in the paragraph 
associated to Eq. (52). The Dirac spinor is defined in Eqs 
( 43|44|45t . 



where x is in fact the center of an overlapping CG cell. 

The evolution of the non-equilibrium distribution po 
is illustrated by a density plot (see Fig. [3| and also by 
surface plot for the final time t = 100 in Fig. [4] We also 
have the data for two further intermediate steps {t = 25 
and t = 75). Fig. [sjand Fig. [6] illustrate, in the same 
manner, the evolution of pi. 



IV. RELAXATION TO QUANTUM 
EQUILIBRIUM FOR A DIRAC PARTICLE IN A 
2D SPHERICAL BOX 

We consider a Dirac particle confined to a spherical 
box of radius R' . Inside the spherical box, the potential 
is negative and constant, outside it is zero: 



C. Discussion 

For that superposition of eight modes with positive 
energy, the relaxation to quantum equilibrium is rather 
fast. In these units, it takes about At = 100 for the 
non-equilibrium distribution to resemble which is 

about one hundred times the Compton wavelength of the 
particle divided by the speed of light. In standard units, 
it is about 10~^^s if we take a Compton wavelength of 
lO-^'^m. 



V{r) 



-\Vo\ r<R' 
r> R' 



(53) 



The corresponding 3D case is analyzed in 



A. Energy eigenstates 

It is similar to the previous case: we can start from 
Eq. (35), replace E hy E ^ \V\ (where |V| is equal to 
1^0 1 oi" equal to zero, depending on whether we consider 
the interior solution or the exterior one) and set = 0. 
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FIG. 4. Surface plot for ip^ip and po at the final time t = 100. 
See Eq. ( |52[) for th e definition of the smooth coarse-graining 
and Eqs ( |43|44|45 ) for the definition of the Dirac spinor. 



Then we have that 

- {{E + \V\)^ - m^)^[ = dl^[ + k{l - k)^[/r^ , (54) 

where k is half-integer (we define n = k — 1/2). 
We distinguish two cases: 

a. = {E ^ \V\Y - rn^ > 0. In that case, the 
equation to solve is 

r^S^^; + r'^n^^lj[ + k{l - k)ilj[ = . (55) 

If we introduce s = vk and = y/s^i^ we find the 
following equation: 



(56) 



which admits Bessel functions of the first and second kind 
as solutions (J^_i(s) and y^_i(s)). 

If we put all the pieces together, we have that 



^P[=VsiaJk-iis)+pY,_.is)) 



(57) 



and 





5 -5 





5 -5 




FIG. 5. Density plots for t G {0,50,100}. The first column 
displays the evolution of the smoothly coarse-grained stan- 
dard density (ip^ip), while the second column displays the evo- 
lution of the smoothly coarse-grained non-equilibrium density 
(pi). The smooth coarse-graining is defined in the paragraph 
associated to Eq. (52). The Dirac spinor is defined in Eqs 
( 43|44|45| ). 





^1 



/^e-^E'e^^^-i^\aJ,_^(s)^f3Y,_^(s)) (58) 



FIG. 6. Surface plot for ip^ip and pi at the final time t — 100. 
See Eq. ( |52|) for th e definition of the smooth coarse-graining 
and Eqs ( |43|44|45 ) for the definition of the Dirac spinor. 
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For tp'2, we have that 



E+\V\+m'2s s 
-in 



.ip[ fc , , 1 , , 1 , , , 

^/x{aJ^_3{s) ^ l3Y^_3{s))) . 



For \Ij2'- 



^2 = -m 



' £;+ IVI +m 



-1 -2fe 



(a4_3(.)+/3r,_3(.))] . (59) 

h. -n^ = {E ^ \V\f - < 0. In that case, the 
equation is 

r^SV; - r^A^V; + k(l - k)ilj[ = . (60) 

If we introduce again s = rn and ijj'^ = v^<I>i, we find 
the fohowing equation: 

s^dl^i + sS.^i - s^^i - n^^i = . (61) 

It admits modified Bessel functions of the first and second 
kind as solutions (denoted by Ik-^ ^k-\)' 
We find that 

4„ = ^e-'^'e^''^-'^)\a'h_.{s) + fi'K^_.{s)) (62) 
and that 



^^2 



E^\V\^m 



1 _ 2k 

. (63) 

The next step is to choose some numerical values for 
the potential Vq and the radius of the spherical box and 
find the energy eigenvalues numerically. In order to do 
that, we impose that the ratio ^ coincides at the bound- 
ary r = R\ 



B. Simulations 

We choose m = 1^ R' = 5 and \Vo\ = 1. We are 
going to superpose energy eigenstates whose eigenvalues 
E satisfy the relations 

< and {E + \V^\f - > . (64) 

Therefore we take 



• Eq. ( [58| ) and Eq. ( 59 ) for the interior solution (and 
we set /3 = for the solutions to be regular at the 
origin) , 



• Eq. (|62| and Eq. (63) for the exterior solution (and 
we set a' = for the solutions to vanish at r = oo.) 

The next step is to find the energy eigenvalues numeri- 
cally, for different values of /c, by imposing that 



V^in,2 



'out,l 



V^out,2 



(65) 



and by looking in the domain defined by Eq. (64). The 
last step is to find p>' by requiring that 



in,l 



'out,l 



(66) 



The guiding spinor used for the simulation is a superposi- 
tion of 6 positive-energy eigenstates with different values 
of k. The details about the spinor can be found in Ap- 
pendix |A] 

For the relaxation simulation, we have considered the 
evolution of the 5 previous non-equilibrium distributions 
for t = 50 and t = 100. The lattice consists of 1024 x 1024 
points distributed uniformly in the box [—3, 3] x [3, 3] and 
the precision on backtracked trajectories is 0.005. If the 
box in divided in 32 x 32 CG cells, the mean percentage 
of good trajectories at t = 100 is 93.22% and the worst 
cell has a percentage of 67.58%. Some of these results 
are illustrated in Fig. ([t]) and Fig. ([8|. 

C. Discussion 

While the relaxation to quantum equilibrium is slower 
in this case (compared to the Dirac oscillator case), it 
must be pointed out that the coarse-graining length is 
smaller here (3/5 of the coarse-graining length used in the 
previous case). Also, if the relaxation time depends on 
the energy spreading, we would expect a longer relaxation 
time in this case anyway. Finally, the Dirac oscillator has 
a strong spin-orbit coupling [22 which may speed up the 
relaxation. 

Another interesting example is the case of a wave- 
function with only positive values of k (which means that 
the direction of rotation doesn't change). In that case, 
non-equilibrium distributions with support around the 
origin will not relax. In order to understand this, we plot 
a few backtracked trajectories (cf Fig. [9| for a superpo- 
sition of the 3 modes with positive values oi k . As we 
can see, trajectories from the inner core are backtracked 
to the inner core, which is a bad thing for relaxation. 
One might argue that only the modes with Jo(x) are 
dominant in that region (these are the eigenstates with 
/c = 1/2 and k = 3/2). We have only tested that the 
inner core is backtracked to the inner core for a superpo- 
sition of 3 modes with positive values of k but it might 
be a general feature of many superpositions with posi- 
tive values of k. For instance, if we simply add one mode 




X 

P3 






FIG. 7. Density plots. The first plot displays the smoothly 
coarse-grained standard density ('0^'^) at t = 100, while the 
5 remaining plots show the smoothly coarse-grained non- 
equilibrium densities at t = 100. The initial non-equilibrium 
distributions (pj, with j G {0,1,2,3,4}) are defined at Eq. 
(|46} and Eq. ( [47| . The guiding spinor -0 is the superposition 
of 6 modes given in Appendix [A] 



with k — —1/2, the core is not disconnected anymore: 
trajectories backtracked from the core can end up in the 
outer region. 



V. CONCLUSION 

We have considered two 2-dimensional systems (the 
Dirac oscillator and the Dirac particle in a spherical 
box) and we have done numerical simulations which show 
an efficient relaxation to quantum equilibrium for Dirac 
fermions, provided that enough modes are superposed in 
order to get sufficient complexity in the dynamics. The 
relaxation takes place despite the absence of nodes thus 
it can be attributed to the intrinsic vorticity of the de 
Broglie-Bohm Dirac velocity field. For the second sys- 
tem, we have also given an example of a non-trivial spinor 
(with positive values of /c) for which there will never be 
complete relaxation. 

This last example does not ruin the idea of relaxation. 
Firstly it is a two-dimensional system. But secondly (and 
more importantly), even in a two-dimensional universe, 
this system does not sit there on its own: it starts to 



FIG. 8. Surface plot for V^t^(t 
(po(t = 0) being defined at Eq. (|46f) 



100) and po(t = 100) 
The guiding spinor 
il) is the superposition of 6 energy eigenstates given in Ap- 
pendix [A] 



interact with other systems, resulting in a more com- 
plex system, where relaxation can take place. However 
this example is interesting with respect to relic non- 
equilibrium particles [27 , which are hypothetical non- 
equilibrium particles from the very early universe that 
would still be in non-equilibrium today. The possibility 
of finding such systems arises because the spatial expan- 
sion in the very early universe stretches non-equilibrium 
length-scales and competes with the natural process of 
relaxation. Then, if these non-equilibrium particles de- 
couple very early and if they are still in non-equilibrium 
at the time of decoupling, they can't be driven to quan- 
tum equilibrium by interaction with other systems. Fur- 
thermore, if they are in an energy eigenstate, the natural 
relaxation does not take place and the non-equilibrium 
can be preserved up to the present day. The last example 
shows that non-equilibrium could also be preserved when 
the spinor is more complex than an energy eigenstate. 
As possible future work, one can study 

• the influence of spatial expansion (which is relevant 
for the early universe and results in a stretching of 
the non-equilibrium length-scales), 

• the role of mass, since it is not fundamental and 
beables should really be attributed to massless 
fermions (discussed in [12 ), 
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tion. For instance, it is still an open question whether 
neutrinos are Dirac or Majorana fermions (a Majorana 
particle being its own anti-particle). The Majorana equa- 
tion is not only relevant for particle physics but also to 
quantum information and condensed matter research (an 
overview of the wide relevance of the Majorana equation 
is given in [28 ). From the point of view of the pilot-wave 
theory, the Majorana equation is also very particular. In- 
deed, in a forthcoming paper [29], we will show that the 
Majorana equation predicts luminal motion for the be- 
able, although the Majorana solutions involve the mass. 
In the same paper, we study relaxation to quantum equi- 
librium for systems governed by the Majorana equation. 
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FIG. 9. Various trajectories originating from different final 
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• relaxation timescales, via the H- function (like in 
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Also the Dirac equation is not the only relativistic equa- 
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Appendix A: Guiding spinors for the second system 

We take |Vb| = 1, m = 1 and i?' = 5. 

For an energy eigenstate, the interior solution is 

i^t{r,e) = «te-*^*e^('=-5)^4_,(K,„r) (Al) 
3 g-iBtgi(fc+i)ei 



E+\Vo\+m 



1 — 2k 

) , (A2) 

riin.T 



where tiin = \/ {E ^ |K)|)^ — m^, while the exterior solu- 
tion (r > R') is 

r^-\r,e) = 4,,e-^^e^(^-^)V'i^,_i (^,.,r) (A3) 

3 ^-iEt^iik+^)e 



m 



-I -2k 



fe_ 1 («o«tr) - p'K^_ s (Koutr)] ( A4) 



The energy eigenvalues are found by matching the ratio 
for the interior and exterior solutions: 

W2 



k 


T7 

rj 


1/2 


0.410077354998218 


3/2 


0.610542082182398 


5/2 


0.812057491976715 


-1/2 


0.598385922365134 


-3/2 


0.356509811273382 


-5/2 


0.510184308650916 



(A5) 



The parameter P' are found by matching for the 
interior and exterior solutions: 



k 




1/2 


-32.6316901377613 


3/2 


-19.79344405979468 


5/2 


-5.13915809445641 


-1/2 


22.59183163168054 


-3/2 


24.1846971959765 


-5/2 


-12.1855792791713 



(A6) 



Finally we take random phases e 



k 




1/2 


0.797881698340871 


3/2 


5.73890975922526 


5/2 


3.97323032474265 


-1/2 


1.74985591686112 


-3/2 


5.11905989575681 


-5/2 


0.61286443954863 



(A7) 



with Kout = Vm^ — E'^. 



The superposition of 6 modes is obtained by summing the 
corresponding eigenstates multiplied by the phase coeffi- 
cients. In the end the spinor is normalized numerically. 
The superposition of 3 (resp. 4) modes is obtained by 
superposing only the first 3 (resp. 4) eigenstates. 



